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Abstract 

The Navier- Stokes equations are solved numerically 
for two-dimensional steady viscous laminar flows. 
The grids axe generated based on the method of De- 
launay triangulation. 

A finite- volume approach is used to discretize the 
conservation law form of the compressible flow equa- 
tions written in terms of primitive variables. A pre- 
conditioning matrix is added to the equations so 
that low Mach number flows can be solved econom- 
ically. The equations are time marched using ei- 
ther an implicit Gauss- Seidel iterative procedure or a 
solver based on a conjugate gradient like method. A 
four color scheme is employed to vectorize the block 
Gauss- Seidel relaxation procedure. This increases the 
memory requirements minimally and decreases the 
computer time spent solving the resulting system of 
equations substantially. A factor of 7.6 speedup in 
the matrix solver is typical for the viscous equations. 

Numerical results axe obtained for inviscid flow 
over a b um p in a channel at subsonic and transonic 
conditions for validation with structured solvers. Vis- 
cous results are computed for developing flow in a 
channel, a symmetric sudden expansion, periodic tan- 
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dem cylinders in a cross-flow, and a four- port valve. 
Comparisons are made with available results obtained 
by other investigators. 


1 Introduction 

The development of a general computer code that 
can predict the flow about complex geometries which 
include complex flow structures is desirable in com- 
putational fluid dynamics. Many numerical schemes 
proposed to date which use the finite difference or 
finite volume formulation of the flow equations were 
written to take advantage of some inherent grid struc- 
ture which then permits flow solutions to be obtained 
efficiently [1] - [3]. A structured mesh can be defined 
as a domain that is discretized such that the neigh- 
borhood of a cell or node can be related to its own in- 
dex number. This structme, which makes the solver 
so efficient, often makes it difficult to obtain reason- 
able grids about complex flow geometries. Many of 
these solution algorithms can be found discussed in 
review papers [4], [5]. 

An unstructured grid flow solver can alleviate 
many of the problems associated with structured 
grids. However, unlik e the structured grid, the cell 
neighborhood of an unstructured grid must be de- 
fined explicitly. The triangle is the simplest and most 
convenient geometric figure that can be used to cover 
a two-dimensional domain. An advantage of using a 
simple triangular shaped cell is the ability to generate 
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grids about arbitrary geometries. Another advantage 
is the ability to add cells in high gradient regions of 
the flow field as well as those regions of the flow that 
are of interest without concern for the surrounding 
cells. The main disadvantages of using an unstruc- 
tured mesh lies in the added complexity and memory 
requirements of the flow solver. 

Several researchers have recently reported success- 
ful results in solving the Euler and Navier-Stokes 
equations on unstructured grids [6] - [11]. The re- 
search reported here considered the use of unstruc- 
tured grids in predicting low Mach number flows 
through internal geometries. To date, there has not 
been much work done toward applying unstructured 
grids to viscous internal flows at low Mach numbers. 

In the research to be described here, the Euler and 
Navier-Stokes equations were discretized on unstruc- 
tured grids composed of triangles in finite volume 
form using primitive variables; however, the conser- 
vation law form was retained. Preconditioning of the 
time derivative term was used to allow efficient cal- 
culations at vanishingly small Mach numbers. The 
equations were solved iteratively using either the im- 
plicit Gauss-Seidel method or an iterative conjugate 
gradient based solver. Several iterative conjugate 
gradient based solvers and matrix preconditioners 
were considered. 

Details of the discretization, the preconditioning, 
the grid generation strategy, and the methods for 
solving the algebraic equations are given in sections 
that follow. Results are given for several test cases 
including the inviscid flow over a bump in a channel 
at subsonic and transonic conditions, viscous devel- 
oping flow in a channel at several Reynolds numbers, 
a symmetric rearward- facing step flow, flow over a 
cascade of cylinders arranged in tandem, and flow 
through a four- port valve. 


2 Governing Equations 

The Navier-Stokes equations were used to model vis- 
cous fluid flow problems in this study. In conserva- 
tion law form and physical coordinates these equa- 
tions can be written in vector form as 


dU_ dE_ <W_ 

dt **" dx dy 


where the vectors 


( 1 ) 


E = 


pu 

pu 2 + P -T xx 
puv — r xy 

(P + e)u — itr rr — vr xy + q x 


and 


F = 


pv 

puv — r xy 

PV 2 + P - Tyy 
(P + e)v — UT X y — VTyy + 


In this work only Newtonian fluids will be considered, 
so the shear stress tensors are defined as 


r xx = + v y) + 2 \iu x - §p(2 u x - v y ), 

T*y = p(Vx +%), ( 2 ) 

Tyy = -\p{u x + V y ) + 2p.V x = Iflfay - U r ), 

and 


Qx — 

q y - -/cT y . 


(3) 


If a further assumption is made that the gas is ideal, 
where p = P/RT with R being the gas constant per 
unit mass, the Navier-Stokes equations can be written 
in nondimensional form as 


with 


dQ(w) dGjw ) dHjw) _ , 

dt dx dy 


P 

u 

V 

f 


0 = 


p 

T 

Pu 

T 

Pv 

T 


L $[(<%- *)r+V + Vi J 


G = 


Pu 

Zf + RP ~ T XX 

Puv 


Z T 


- T a 


xy 


L T 


f (CpT + % + X )U - UT XX - VT xy - CT X J 


and 





Pv 


p 


. T 
Puv ~ 

u = 

pu 

pv 

H = 

~f T *y 

rf + RP-r yy 


e 


^(CpT + Y + y)^ “ V>T X y — VTyy ~~ CTy 


2 


where 


= |^f ( 2u * ~ v v)> 

= f(*x + V>, 

T yy =.§7|f( 2 % _u x), 

C = 


T XX 


T *y 


(5) 


CpjxFt 
Pr Re ' 


The quantity M is the reference Mach number, de- 
fined as Uref/y/lRTrtf, which appears when the 
equations are cast in nondimensional form. Here and 
in the matrices that follow V is used for convenience 
to represents the quantity 7 — 1 . This equation can 
be rewritten as 


The Reynolds and Prandtl numbers are defined as 


Pre/^re/^rt/ 

f*ref 


Pr = 
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where the Jacobian matrices, A ( and A x can be writ- 
ten as 


respectively. The refers to a term that is nondi- 
mensional. The variables subscripted ref are refer- 
ence quantities specific to a particular flow calcula- 
tion. The test cases presented in this work involve the 
laminar flow of air where the viscosity is assumed to 
follow the Sutherland formula. It is important to note 
that although these equations are written in terms of 
primitive variables(P, u,v,T), they are still in conser- 
vation law form. 

All subsequent equations are nondimensional so the 
" is dropped for convenience. 

3 P reco ndit io ning 

Solving the compressible flow equations for low Mach 
number cases is difficult because the resulting system 
is stiff due to the large ratio of acoustic to convective 
velocities. A temporal preconditioning is used in this 
work to remove this stiffness. 

The approach will be demonstrated by using the 
Navier-Stokes equations in one dimension. The 
nondimensional form of these equations is written as 
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( 6 ) 


It is important to note that the quantity P/RT is 
well behaved as Mach number goes to zero since it 
is equal to the nondimensional density, pj p r ef • Also 
note that R = (7 M 2 )"* 1 . 

The first column of the Jacobian matrices contains 
the coefficients associated with pressure. In A t , these 
coefficients go to zero as the Mach number goes to 
zero. This results in the acoustic time scale restric- 
tion associated with pressure. As the Mach num- 
ber approaches zero, a vanishingly small time step 
is needed to keep the coefficients multiplying pres- 
sure finite. More importantly, the system of equa- 
tions become singular as the Mach number goes to 
zero, as will be discussed below. For a finite time 
step, the time derivative of pressure will vanish from 
the equations. In the limit as the Mach number ap- 
proaches zero, the equations reduce to their incom- 
pressible form. Since the time derivative of pressure 
does not appear in this form of the incompressible 
equations, the equations contain no pressure history. 
There is a more mathematical way of looking at this 
problem. The system can be rewritten in the form 

dw . . dw . _! dG v {w ) 

ar +A ' A 'to= A ' -fc- (8) 
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As Mach number, Af, becomes small, A t becomes 
ill-conditioned, i.e., the determinant of At becomes 
small and errors due to round off error will become 
large when computing At” 1 . In the limit as M goes 
to zero, At -1 is unbounded and the system is singu- 
lar. The result is often slow convergence due to this 
stiffness when using a compressible code to compute 
a flow at very low Mach number. In addition, the 
eigenvalues(f/ + C, U — C, and I7)of the Jacobian ma- 
trix become farther apart as the Mach number goes 
to zero. A remedy for this is write the equations in 
the form 


A p 


dw . dw dw 

d^ + A, 'm +At d^- 


dGj w) 

dx 


( 9 ) 


The preconditioning Jacobian matrix, A p , is defined 
as 


f i. 0 1 

A 5 I _E 

Ap — T PT HT 7 

1 , M 3 Tu 3 y M*TPu yM*rPu_ 

y “T 2T T 2T 5 

This Jacobian matrix used in the preconditioning is 
of the same form as the matrix At, but the depen- 
dence of Mach number is removed from the terms 
that are causing the ill-conditioning. This essentially 
attempts to cluster the eigenvalues around the con- 
vective speed. It may be possible to simplify Ap by 
setting some of the nondiagonal terms in columns two 
and three equal to zero. 

For a time dependent calculation at low Mach num- 
ber, the preconditioned equations are advanced in the 
pseudo time frame as well as the real time frame. 
At each physical time step, the equations are it- 
erated to convergence in pseudo time, r. At this 
point the pseudo time term vanishes and the time 
dependent Navier-Stokes equations are satisfied. The 
pseudo time iterations also remove the linearization 
error from the solution at each physical time level. 
If the low Mach number flow computation is steady, 
it is only necessary to integrate the equations in the 
pseudo time frame. This is done by setting the phys- 
ical time step to a very large number to remove the 
effect of the physical time derivative from the precon- 
ditioned equations. 

At higher Mach numbers the pseudo time term is 
not needed, although convergence does not appear to 
deteriorate with its continued use. 

Pseudo time terms for the full two-dimensional 
equations are added to the diagonal blocks of the 
sparse matrix, A, in Eq. (18). These terms are 
formed by a direct extension of the one-dimensional 
example above. 


4 Grid Generation 

The triangular shaped cell is the simplest geometric 
shape that can be used to cover a two-dimensional 
computational domain. With an unstructured grid, 
individual cells can no longer refer to their neighbors 
simply by incrementing an index as in a structured 
grid. Instead, the neighborhood of a cell is deter- 
mined through a connectivity matrix. This connec- 
tivity matrix usually contains cell based information 
as well as edge based information. Details of the con- 
nectivity matrix required by the computer flow code 
developed in this work will be discussed later in this 
section. 

The use of a triangular unstructured grid formula- 
tion has some distinct advantages over a structured 
grid. One advantage to using triangular cells is that 
with them it is easy to generate grids about complex 
geometries. This reduces the amount of time required 
to generate a suitable grid. Also grid adaptation can 
be done locally without adding unnecessary cells to 
other regions of the domain. 

The method used in this work follows a path sim- 
ilar to that of Holmes and Snyder [12] to triangu- 
late a region. First, the boundaries that describe the 
computational domain are defined. Next, the bound- 
aries are discretized in a counter-clockwise direction. 
These discrete points axe then triangulated using the 
Delaunay criterion that no other point in the com- 
putational domain lies within the circumcircle of a 
given triangular cell. Points must now be added to 
the domain to obtain a reasonable grid. A new point 
can be added based on any criteria one chooses. The 
grid point insertion in the current work is done ac- 
cording to one of the following three geometric cri- 
teria: improve the triangle with the smallest aspect 
ratio, reduce the maximum area triangle, reduce the 
size of the triangle with the largest circumcircle ra- 
dius. These geometric constraints can be used in any 
combination. Some other criteria that could be used 
for local retriangulation are increase minimum angle, 
decrease maximum angle, and maintain equal length 
sides, to name a few. An example of a coarse grid 
about a square hole generated in the above manner 
is shown in Fig. 1. 

The flow code requirements dictate the type of out- 
put that the grid generation scheme must provide. A 
connectivity array must be generated for an unstruc- 
tured grid so that a cell neighborhood is completely 
defined for the flow code. The code in the current 
work was based on a cell centered scheme. Here the 
triangle itself was the control volume used in the fi- 
nite volume formulation. 

Connectivity is determined by cell nodes, cell faces, 
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Figure 1: Final coarse triangulation of domain 


and face cells. Cell nodes are the nodes at the ver- 
tices of a triangle. The cell faces are the edges of 
a triangle. The face cells are the cells that lie on 
either side of a given edge. Typically the x, y coordi- 
nates of the cell nodes are the only grid floating point 
numbers required as input by a flow code. The grid 
connectivity is defined through integer arrays. Two 
connectivity arrays were needed in this work to define 
the discretized geometry of a computational domain, 
Fig. 2. The two-dimensional array NCELL contains 
the edge and node connectivity for a given cell. The 
first dimension of NCELL contains 6 elements. The 
second dimension has a length n, where n is the to- 
tal number of triangular cells in the computational 
domain. Consider the cell number i. The first three 
elements of the first dimension of array N CELL are 
the edge numbers of cell i. The last three elements 
of the first dimension of array NCELL are the ver- 
tex node numbers of cell i. This allows the access 
to the edge and node numbers that define a given 
cell, in this case cell number i . The array N FACE 
is also two-dimensional. It contains cell connectivity 
for a specific edge. The first dimension is of length 
2. The second dimension is of length to, where to is 
the total number of edges that make up the compu- 
tational domain. The second dimension identifies the 
edge(in this example, j). The first dimension of array 
N FACE contains the cell numbers that are adjacent 
to one another sharing the common edge j. These two 
integer arrays along with the floating point arrays x 
and y define the geometry for the present scheme. 

Boundary information must be defined explicitly. 
Solid wall, exit, and inlet boundaries are implied 



Figure 2: Connectivity requirements for a single cell 


through the edge connectivity array, NF ACE. For a 
solid wall boundary, one of the elements of the first di- 
mension of the array NF ACE will contain the value 
0. This tells any cell that refers to that edge that 
it borders a solid wall boundary. Sunil arly, an exit 
boundary is adjacent to a cell number of —1, and an 
inlet boundary borders a cell number of —2. Periodic 
and symmetric boundaries are handled through spe- 
cial connectivity. This again is done by including the 
appropriate cell information in array N FACE. A 
symmetry boundary cell will have an edge that bor- 
ders itself. So the first dimension of array N F ACE 
for the symmetry face will have both elements refer- 
ring to the same cell number. For a periodic bound- 
ary the elements will refer to cell numbers that are 
separated by one periodic pitch. Boundary conditions 
will be addressed later. 


5 Discretization Technique 

The finite volume formulation of the governing equa- 
tions is well suited for application to an unstructured 
discretization of the flow domain. The nondimen- 
sional Navier-Stokes equation written in differential 
form Eq. (4) is first recast in integral form for an 
arbitrary volume, V as 



dV = 0. 


( 10 ) 


Using Gauss’s theorem, the area integral of the flux 
derivatives can be rewritten as the surface integral of 
the flux quantities around the perimeter of the vol- 
ume, V. This allows the Eq. (10) to be written as 


d_ 

dt 


J QdV + j ( Gdy - Hdx) 


= 0 . 


( 11 ) 
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For each control volume consisting of a triangular el- 
ement, Eq. (11) is evaluated as 

TMiQi) + £(<?, A Vj - A Xj) = 0. (12) 

dt U 

where Qi is the vector of conserved quantities in cell 

i, Gj and Hj are the flux vector quantities across edge 

j, and Axj and Ayy are the differences in Cartesian 
nodal coordinates that define edge j . The summation 
on j proceeds in a counterclockwise manner around 
the edges of cell i. Also it is understood that Axj 
stands for x{end) — x (beginning) as the evaluation 
proceeds in a counterclockwise manner around the 
sides of a control volume. The quantity Ai is the 
area of cell i defined as 


Ai = y/s($ — a)($ — 6)(s - c). (13) 

The variable s is the semiperimeter of cell i , 

s =i(a + 6 + c) (14) 

and the quantities a, 6, c refer to the lengths of the 
sides of the cell i. Cell face flow quantities required 
by Eq. (12) for the computation of the inviscid flux 
terms were approximated by using the average of the 
cell centered values on both sides of a given cell face. 
The numerical integration of these quantities around 
the edges of the cell results in a central difference 
scheme that is second order accurate in space. The 
viscous terms required the computation of the deriva- 
tives on the faces of the triangle control volume. To 
compute these terms, the level 2 cells shown in Fig. 3 
were used and a different path integral was evaluated. 
Again this yields a second order accurate scheme in 
space. A total of 10 cell centered quantities was used 
in the computation of the viscous quantities of the 
summation term of Eq. (12). Specifically the viscous 
terms in the x- momentum equation are written as 

3 

^(-r xx Aj/+t, s A x ) ; - 
j= i 



where the summation is over the three sides of the 
triangular control volume. A typical derivative can 
be recast in integral form as 



Figure 3: Cell level dependence 



where S' is the sum the areas of the two cells across 
a given edge, and the integral is along the path that 
traverses the outer boundary of the two cells in a 
counterclockwise direction. These derivatives are in- 
terpreted as mean values over the area S ' . 

The system of equations was integrated in tune us- 
ing an implicit scheme written in delta form. Newton 
linearization was used on nonlinear terms. For exam- 
ple, the terms 


P 

T 

Pu 

T 

Pv 

T 


PI P 

~4 + 4ap-4-at, 

T T p2 

~A + iA« + ^-AP-£^AT, 

p p p p2 

Pv P . v . _ Pv 

~ — + -r Av + -r A P — -AT 

p p p p2 


(17) 


are substituted into the continuity and momentum 
equations. The ~ terms take on provisional values of 
the primitive variables; and the delta quantities rep- 
resent the differences between values at the new time 
level and the provisional values, e.g., A u = u — u. 
These equations can be iterated at a specific point 
in time until the linearization errors are reduced to 
a satisfactory level. For steady state problems, this 
iteration process at each time level is not necessary. 
Similar terms are used to linearize the energy equa- 
tion. However, the nonlinear dissipation terms in the 
energy equation were linearized by evaluating them 
explicitly in terms of the provisional values, rather 
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than using Newton linearization. The result of the 
above substitutions is the matrix equation 

Ax — b. (18) 

The matrix A contains the linearized terms which 
multiply the vector of unknown delta quantities, x. 
The vector b represents the residual of the equations 
which should go to zero as the solution approaches 
convergence. There are four equations written for 
each cell in the computational domain. For the cen- 
tral difference formulation given here, each cell is de- 
pendent on the level one cells through the convective 
terms as well as the level two cells through its diffu- 
sive terms, as illustrated in Fig. 3. This along with 
the unstructured grid yields a sparse matrix whose 
elements are block 4x4 matrices. The general row of 
the matrix A has ten non-zero blocks of coefficients. 
The number of blocks vary near a boundary. 


6 Artificial Dissipation 

Artificial dissipation was needed in the current imple- 
mentation of the flow equations to prevent the odd- 
even decoupling seen in central difference computer 
flow codes. Two schemes were used in this work. 
The first scheme was based on the research of Jame- 
son and Mavriplis [6] and was used specifically for 
inviscid subsonic and transonic test cases. The sec- 
ond version was developed to be used with the low 
Mach number test cases where preconditioning was 
employed. Both dissipation schemes were added ex- 
plicitly to the system of flow equations. 

The second type of dissipation was developed to 
be used with the temporally preconditioned scheme. 
This form of the dissipation was used for the vis- 
cous subsonic flows computed in the present work. 
Here a biharmonic operator was used on the primitive 
flow variables. The second difference of the primitives 
were computed as 


3 

a Pi2 = Y. [Wj] - 3u;i, 

J= 1 

where w was the vector of primitive variables. Again 
the summation was done over the index j. The result 
represents the second difference of the variables in 
cell i. The fourth difference is then computed by 
summing the third difference over the edges of the 
cell, 


3 

a PiA = Y ” 3a P*2- 

J = 1 


The resulting fourth difference was then premulti- 
plied by the preconditioning matrix A p . Then it was 
multiplied by the appropriate coefficient to make it 
consistent with the other terms and included explic- 
itly on the right-hand-side of the system of equations. 


7 Boundary Conditions 


For the subsonic viscous flow cases in this study, the 
inlet boundary conditions were imposed by specifying 
the velocity components, u and u, as well as the static 
temperature, T. The static pressure was extrapolated 
from the interior to the inlet. At the exit, u, v, and T 
were extrapolated downstream. The static pressure, 
Z 5 , was specified. At the solid wall the static tem- 
perature was specified. The u and v velocity com- 
ponents were set to zero to enforce the no-slip con- 
dition. Static pressure was specified as symmetric 
in the ghost cell to give a zero normal derivative at 
the solid surface. Symmetry and periodic boundary 
conditions were imposed by simply specifying the ap- 
propriate cell connectivity. At a symmetry boundary 
cell, values were reflected across the boundary. At a 
periodic boundary cell, values were transposed by the 
periodic pitch of the computational domain. 

For supersonic viscous flow, boundary conditions 
at a solid wall remained the same as in subsonic flow. 
At the inlet, all flow quantities, P, u y v y and T were 
specified. At the exit all the flow quantities were 
extrapolated. 

It is important to note that one of the test cases 
computed to verify the code was inviscid. For this, 
the viscous boundary conditions described above 
were modified appropriately. Since the inflow and 
outflow were subsonic, it was only necessary to mod- 
ify the no-slip solid wall boundary condition to a tan- 
gency boundary condition. 

In the current work, temporal preconditioning was 
used to compute low Mach number flows as described 
in the previous section. The preconditioning did not 
affect how the boundary conditions were specified, 
but it is interesting to note that the eigenvalues of 
the system were modified. The eigenvalues must now 
be obtained from the matrix that results from the 
product A p “ 1 A r . In one dimension 


A p 


-l 


—jM 2 Tu 7 

uRT RT n 

___ p u 

{Tti a -2RT)T TTu yPT 

— 2 p r P J 


where T = 7 — 1 . One eigenvalue remains unchanged 
that is Ai = U . The other two eigenvalues A2,3 for 
the given one-dimensional system take on a similar 
but yet a more complex form than that described 
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by Withington et al. [15]. However, in the limit as 
Mach number approaches zero, the eigenvalues of the 
system are 


Ai = U 

and 

. U ± y/U 2 + 4 7 T 
A 2 ,s= 2 * 

By substituting in the appropriate nondimensional 
quantities the resulting ratio of largest to smallest 
eigenvalue takes on the value of (1 4- v / 5)/(l — y/E). 
This is the same quantity that was obtained by With- 
ington et al. [15] with the ratio of specific heats set 
equal to one in the present work. The preconditioning 
essentially allows all of the equations of the system to 
be integrated at the same pseudo- time rate. This can 
be compared with a scheme without preconditioning 
where the ratio of largest to smallest eigenvalues is 
infinity. 

8 Sparse Matrix Solvers 

The system of algebraic equations being solved in 
the present implicit unstructured grid formulation is 
represented by Eq. (18). The matrix A is sparse. 
There is usually no particular pattern to the nonzero 
elements when the matrix arises from an unstruc- 
tured grid formulation. However, the blocks on the 
main diagonal of this sparse matrix always have some 
nonzero entries. Figure 4 shows a representative form 
of the sparse matrix A. The solid squares repre- 
sent 4x4 blocks with at least the diagonal elements 
of the diagonal block being nonzero. The remain- 
ing blocks contain zeros. In the present method, a 
two-dimensional viscous flow computation requires a 
maximum of ten 4x4 blocks in each row of the A 
matrix. In contrast, an implicit structured solver can 
be written such that the resulting A matrix on the 
left-hand side has some special structure that allows 
the matrix equation to be solved by some well es- 
tablished methods. The equations can often be cast 
in a form that results in a block bidiagonal or block 
tridiagonal matrix. This structure is not generally 
available to the solution of the flow equations written 
for am unstructured grid. 

Several iterative methods were examined in this 
work. The first was a point Gauss-Seidel scheme 
where only the diagonal elements of the diagonal 
blocks of matrix A were retained on the left-hand side 
as unknowns. This scheme was successful for many 
of the simpler problems but was prone to divergence 
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Figure 4: Form of sparse matrix A 


when starting with poor initial conditions. It seemed 
to be very sensitive to lack of diagonal dominance. 

Another iterative scheme used was the point block 
Gauss-Seidel method. Here the diagonal 4x4 blocks 
of matrix A were retained on the left hand side. The 
remaining matrix equation was solved using LU de- 
composition. The L and U matrices refer to the lower 
and upper triangular decomposition of the diagonal 
block of the matrix A. This was found to be more 
robust than the previous scheme. 

Even though the full sparse block matrix is N x 
JV, it is only necessary to store the nonzero blocks. 
This gives a maximum block matrix of 10 x N for a 
viscous code. However, the bandwidth could still be 
the maximum, N. 

The commercially available sparse iterative solver, 
SITRSOL [17], which resides on the Cray YMP as 
a callable subroutine was also evaluated for solving 
the above matrix equation. SITRSOL takes advan- 
tage of the matrix sparseness by only storing the 
nonzero entries. The package makes available to the 
user several iterative methods as well as precondi- 
tioners for solving non-symmetric positive indefinite 
sparse linear systems. In the present work the bi- 
conjugate gradient method, the generalized minimal 
residual method, and the generalized conjugate resid- 
ual method were considered. An incomplete LU 
preconditioner was also used. These three iterative 
methods are of the preconditioned conjugate gradi- 
ent type. 

The iterative solver SITRSOL was used on one of 
the test cases to be shown in the results section and 
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its effectiveness was compared with that of the point 
block Gauss-Seidel method. The conclusions shown 
in this work are provisional. More experience needs to 
be obtained to make a true evaluation of the various 
solvers and preconditioners. 

The solution of Eq. (18) using a point block Gauss- 
Seidel method suffers from recurrence. The penalty 
is seen in vectorization. This recurrence can be elim- 
inated with a minimum effect on the solution con- 
vergence rate by using a coloring scheme. The idea 
comes from a problem which arose in graph theory. A 
theorem states that a map can be colored with only 
four colors such that no two regions of the same color 
share a border. The conjecture was proven through 
exhaustive computation by Appel and Haken [16] in 
1976. This theorem was implemented by first color- 
ing the unstructured grid according to the theorem 
and storing all cell numbers of given color in an in- 
teger array. The scheme was most efficient when the 
number of cells in each color integer array was about 
equal. 

This gave a color map that was then used as input 
to the flow code. The Gauss-Seidel algorithm was 
then written to contain four loops corresponding to 
the four colors of the colored grid. Each single col- 
ored loop contained no level 1 cell recurrence, so it 
was vectorized. On a typical problem in the present 
study, the solution time for the algebraic system(the 
Gauss-Seidel subroutine) was reduced by a factor of 
7.6 times by using this four color partitioning. Re- 
currence is still present but only through the level 
2 cells, illustrated in Fig. 3, required in the viscous 
terms. The result is that the quantities in the level 2 
cells are lagged from the previous iteration time step. 
However, this does not seem to effect the convergence 
rate. 

9 Results 

The results presented in this section will be used to 
demonstrate two conclusions. The initial results will 
show the validity of the code. And later results will 
indicate the versatility of the unstructured grid over 
the structured grid formulation. Comparisons will be 
made with data available from other investigators. 

9.1 Bump on Wall 

Inviscid flow over a bump in a channel was computed 
at two values of Mach number. A Mach number of 
0.5 was used for the first test case. Figure 5 shows 
Mach number contours. The flow was subsonic so 
the inviscid flow was symmetric about the middle of 
the bump. This could be seen more clearly when the 



Figure 5: Constant Mach number contours for flow 
over a symmetrical bump in a channel, M xn = 0.5 



0.5 1.0 1.5 2.0 2.5 


Figure 6: Upper and lower wall Mach number distri- 
bution 

upper and lower wall Mach number distributions were 
plotted. The results of this subsonic case compared 
well with those reported by Ni, [18]. 

A second test case was computed at a Mach num- 
ber of 0.675 at the inlet. Here the flow was transonic 
over the bump. The grid used for this test case can be 
seen in Fig. 7. A supersonic bubble was formed on 
the bump, Fig. 8* The location of the shock was 
shown clearly in the plot of upper and lower wail 
Mach number distributions, Fig. 9. The location of 
the shock compared well with the results of Ni, [18] 
as well as that of Chima et al, [19]. The sonic line 
that impinged on the aft side of the bump w r as at a 
distance of 72 percent of the chord length from the 
head of the bump in the above cases. The present 
case locates the sonic line at 73 percent of the chord 
length. 

These inviscid test cases required the addition of 
artificial dissipation. For the subsonic case, a fourth 
difference was added to prevent the odd-even de- 
coupling of the solution seen in central difference 
schemes. The transonic case also required the addi- 
tional second difference to prevent oscillations from 
occurring about the discontinuity. In both cases the 
dissipation model was similar to that of Jameson et 

al. [3]. Later Jameson and Mavriplis [6] implemented 

this type of dissipation model for an explicit unstruc- 
tured grid flow solver. 
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9.2 Developing Channel Flow 



Figure 7: Computational grid for the symmetrical 
bump in a channel test case, Af,* n = 0.675 



Figure 8: Constant Mach number contours for flow 
over a symmetrical bump in a channel, Mi n = 0.675 
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Figure 9: Upper and lower wall Mach number distri- 
bution 



Developing flow in a channel was used to validate the 
code for viscous flows. It also served the purpose of 
testing the preconditioning used for computing low 
Mach number flows. Comparisons were made be- 
tween the Gauss-Seidel method and the solver SITR- 
SOL for solving the sparse matrix equation. 

The code was validated on four developing channel 
flow test cases. A low inlet Mach number flow of 0.05 
was used to compute flows at Reynolds numbers of 
1, 20, 150, and 1500 based on the inlet uniform ve- 
locity, density, and full channel height. Because the 
inlet Mach number was held constant, the channel 
height was varied to obtain the appropriate Reynolds 
number. Unstructured grids of 1114,1969,4800, and 
4800 cells were used for the Reynolds number flows 
of 1, 20, 150, and 1500 respectively. Uniform flow 
enters the channel with a nondimensional uniform 
velocity of one and accelerates to a nondimensional 
centerline velocity of 1.5. In order to compare with 
published results for incompressible flows, it was nec- 
essary to make a correction to the centerline velocity 
at the low Reynolds numbers due to the larger den- 
sity variation from the inlet to the exit of the chan- 
nel. Figure 10 shows the centerline velocity of the 
channel flow at various Reynolds numbers. These 
results were compared with other computations by 
Tenpas and Pletcher [20], Morihara and Cheng [21], 
and Chilukuri and Pletcher [22]. At a Reynolds num- 
ber of 20 the centerline velocity of the current study 
slightly under predicted the centerline velocities ob- 
tained by the other investigators near the exit of the 
channel. At a Reynolds number of 1500, the results 
of the present study show a more rapid acceleration 
of the flow than indicated by the solution of the par- 
tially parabolized Navier- Stokes equations. 

Typical convergence histories for the code are 
shown in Fig. 11. The convergence criteria was based 
on the residual of the continuity equation in delta 
form which should approach machine zero as the solu- 
tion goes to a steady state. The solution of the matrix 
equation was done by the block Gauss-Seidel method. 
In general, the solution converged at nearly the same 
rate over a wide range of Mach numbers holding the 
Reynolds number equal to 20 for the four flow test 
cases. This illustrates the benefits of the precondi- 
tioning. Without preconditioning, it was necessary 
to run the code at a much smaller time step thus de- 
creasing the rate of convergence. At Mach numbers 
lower than 0.1 the code without preconditioning did 
not converge. 

The sparse matrix solver SITRSOL was used for 
comparison with the Gauss-Seidel method. The con- 
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Figure 10: Centerline velocity profiles for develop- 
ing flow in a channel at M in = 0.05 with Rt'h = 
1,20, 150,1500 



Figure 11: Convergence history for developing chan- 
nel flow over a range of Mach numbers at Rch = 20 
using the block Gauss-Seidel solver 



Iteration Count 


Figure 12: Convergence history of developing flow in 
a channel, Fith = 2(3 computed with sparse matrix 
solver 

vergence history for three different conjugate gradi- 
ent like methods with the ILU used as a precondi- 
tioner is shown in Fig. 12. The Gauss-Seidel method 
without the coloring algorithm took 13.5 minutes on 
the Cray YMP. The same computation with a color 
map supplied for vectorization of the Gauss-Seidel 
algorithm gave a speedup factor of 7 .6 over the stan- 
dard Gauss-Seidel matrix solver. The overall com- 
puter time was reduced to 11.4 minutes, or a speed 
up of 15.5 percent. This suggests that more attention 
should be given to the vectorization of other parts 
of the flow code. The coloring scheme did not have 
much effect on the convergence history of this viscous 
calculation. The same grid was used to make com- 
parisons with SITRSOL. The bi-conjugate gradient 
method took 9.3 minutes of computer time to reach 
about the same level of convergence as the Gauss- 
Seidel method. The generalized conjugate gradient 
residual method required 10.23 minutes of computer 
time. About the same level of convergence was ob- 
tained by the generalized minimum residual method 
in 5.5 minutes. 

9.3 Sudden Expansion 

The previous test cases could have easily^been com- 
puted using a structured grid approach. The sud- 
den expansion test case demonstrates the capabil- 
ity of the unstructured grid generation and its abil- 
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ity to obtain a grid in a domain that would other- 
wise need a patched or masked grid to work for a 
structured flow code. The results from this com- 
putation were compared with the experimental re- 
sults obtained by Durst et al. [23]. They noted that 
though at lower Reynolds numbers the flow was sym- 
metric about the centerline of the expansion, there 
were three-dimensional effects near the separated re- 
gions. A plane symmetric sudden expansion with a 
downstream channel height to step height ratio of 3:1 
was computed. The Reynolds number for this flow 
was 56 based on the upstream channel height and 
the centerline upstream velocity. A fully developed 
parabolic profile was prescribed at the inlet which was 
located one step height upstream of the expansion. 
The Reynolds number was computed at a streamwise 
location 0.25 step heights upstream of the expansion. 

It was interesting to note that the profile at this lo- 
cation was already anticipating the expansion corner. 
The flow neax the wall begins to accelerate; and to 
conserve mass, the centerline velocity decreases. The 
flow separates at the step, reattaches downstream, 
and returns to a fully developed profile about ten step 
heights from the expansion. The streamwise velocity 
component at six specific channel locations are shown 
in Fig. 13. Comparisons were made with the laser 
anemometer experimental data presented by Durst et 
al. [23]. The centerline velocity distribution was com- 
pared with the laser anemometer data and with the 
viscous-inviscid interaction computational method of 
Kwon et al. [24] and is shown in Fig. 14. The pre- 
dicted centerline velocity appears larger than the ex- 
perimental values downstream, but the correct value 
of one-third the upstream fully developed centerline 
velocity was obtained in the present calculation. 

9,4 Periodic Tandem Circular Cylin- 
ders in Cross Flow 

The flow was computed over a cascade of tandem 15 
circular cylinders. This computation should be of 
practiced interest in that geometries of this sort are 
encountered when modeling flow through heat ex- 
changers. These tube heat exchangers can be found 
in automobile radiators, room heaters and gas and 
air heaters. With the unstructured grid formulation, 
it was easier to generate a computational grid about 
in-line as well as staggered cascades of tubes. Some of 
the geometric quantities that affect the flow charac- 
teristics of the heat exchanger are the size and shape 
of the tubes as well as their vertical and horizontal 
spacing. This type of parametric study is ideal for 
the unstructured grid formulation. 

The model problems presented here were compared 



Figure 13: Velocity profiles for a laminar flow in 
a channel with a 3:1 symmetric sudden expansion, 
Ren = 56 
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Figure 14: Centerline velocity distribution for a lam- 
inar flow in a channel with a 3:1 symmetric sudden 
expansion, Re & = 56 
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Figure 16: Periodic tandem circular cylinders in cross 
flow, Re r — 1 

with the incompressible numerical results of Gordon 
[25]. Uniform flow conditions were prescribed per- 
pendicular to the cascade upstream of the first bank 
of tubes. Periodic boundary conditions were imposed 
at the upper and lower geometric boundaries to simu- 
late an infinite number of parallel rows. The tubes in 
this case were in-line. An upstream Mach of 0.05 w T as 
used for both computed test cases. Both cases were 
computed on the same geometry. The geometry used 
in both test cases and the grid used for the second 
test case is shown in Fig. 15. 

The first test case was computed at a Reynolds 
number of 1.0 based on the upstream conditions and 
cylinder ra dius . The velocity vectors are shown in 
Fig. 16 and compare well with the streamlines com- 
puted by Gordon [25]. At this Reynolds number the 
flow was nearly symmetric about both cylinders indi- 
cating that there was minimal influence of one bank 
of cylinders on the other. 

The second test case was computed at a Reynolds 
number of 20.0 based on the same conditions as the 
test case above. Here the flow separates behind both 
cylinder banks. An enlargement of the velocity vec- 
tors near the cylinders is shown in Fig. 17 and the 
length of the separated regions behind both of the 
cylinders compares well with those computed by Gor- 
don [25]. As noted by Gordon [25], the length of the 
separation behind the second bank of tubes is slightly 
smaller than that behind the first bank. The first set 
of cylinders accelerates the flow in the freestream so 


Figure 17: Periodic tandem circular cylinders in cross 
flow, Re r = 20 

the slower wake flow impacts the second set of cylin- 
ders. The flow was symmetric about an imaginary 
horizontal line that passed through the centers of the 
cylinders. There does not seem to be any difference 
in the angular location of the actual separation point 
on either cylinder. 

9,5 Four Port Valve 

The final results are presented to show the versatil- 
ity of applying boundary conditions when using an 
unstructured grid flow solver. The geometry repre- 
sents a two-dimensional valve with four ports where 
inlet or outlet flow boundary conditions can be speci- 
fied. A reference was made to such a flow geometry in 
an article by Ackert [26]. The actual flow conditions 
were not given. Here the flow was computed at two 
Reynolds numbers. Fluid enters through a channel on 
the left, enters the circular cavity, and exits through 
a channel at the bottom. For both test cases fully de- 
veloped conditions were prescribed at the inlet. The 
flow' redevelops along the open channel. Both test 
cases were computed with an inlet Mach number of 
0.05. An interesting aspect of the geometry was that 
the closed valve ports acted as driven cavities. The 
unstructured grid formulation allows the application 
of exit boundary conditions at any or all of the three 
remaining ports. This type of valve geometry can be 
found in an application like fluid networks. 

Fluid flow was first simulated in the valve geome- 
try at a Reynolds number of 10 based on inlet con- 
ditions and channel height. The velocity vectors of 
this steady state flow are shown in Fig. 18. Here the 
fluid near the wall of the inlet channel was acceler- 
ated as the corner of the cavity w^as anticipated. A 
clockwise rotation of the fluid was followed through 
the circular volume. The fluid in the closed cavity 
ports was driven in a counterclockwise rotation. The 
band of fluid then enters the open lower exit channel 
and again becomes fully developed. 
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Figure 19: Four port valve, Ren = 50 

Next the same valve geometry was used to simulate 
the fluid flow at a Reynolds number of 50 based on 
inlet conditions and channel height. The velocity vec- 
tors of this computation are shown in Fig. 19. Again 
a fully developed velocity profile was specified at the 
inlet to the channel. The velocity of the fluid near the 
wall accelerates as it approaches the entrance to the 
circular chamber. Contrary to the previous case, the 
banded fluid actually drove a large volume of fluid in 
a counterclockwise direction. This had the effect of 
driving the closed valve ports in an opposite rotation 
direction from that of the lower Reynolds number test 
case. Also, the band of fluid did not diffuse as much 
across the circular volume. In the open exit channel 
the fluid redevelops into a fully developed parabolic 
profile. 


10 Conclusions 

A two-dimensional unstructured grid implicit flow 
solver was described. Although only internal flow 
problems were considered in this study, the method 
is believed to be applicable to external flows as well. 
The compressible flow equations were discretized in 
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finite volume form. 

Several conclusions can be drawn from the present 
study. 

1. A triangular unstructured grid can be generated 
about very complex geometries where the use of 
a single structured grid cannot be considered in 
most cases. This gives the advantage that a sin- 
gle computer code can be used in a wide vari- 
ety of flow geometry applications. However, this 
advantage is somewhat dampened by the com- 
plexity of coding required for solving a system of 
differential equations on this unstructured grid. 
Details such as boundary conditions are more 
difficult to implement on an unstructured grid. 

2. It was found that the diagonal block Gauss- 
Seidel solver was more robust than the point di- 
agonal Gauss-Seidel version of solver. The diag- 
onal point solver seemed sensitive to initial con- 
ditions and diagonal dominance. 

3. A coloring scheme was used to take advantage 
of the vectorization of the implicit Gauss-Seidel 
solver. A minimum of extra storage was neces- 
sary for a significant reduction in computer time. 
The time spent in the solver was decreased by a 
factor of 7.6. It was found that the recurrence in 
the viscous fluxes had little affect on the conver- 
gence of the solution to a steady state. 

4. The use of the sparse matrix iterative solver al- 
lowed a much larger time step to be used than 
that of the Gauss-Seidel solver. However, every 
time step using the sparse matrix solver was sig- 
nificantly more expensive. Even so, the sparse 
solver ran at 2 to 2.5 times faster than the block 
Gauss-Seidel solver. Several different conjugate 
gradient like solvers were tested with the matrix 
preconditioners available with SITRSOL. Some 
of the preconditioners did not allow a solution 
to the equations. The incomplete LU precondi- 
tioner was found to be the best. The solvers all 
exhibited the same basic convergence rate. The 
difference in the solvers was in the time that was 
required to obtain the same convergence level. 
The generalized minimum residual method was 
found to be the fastest for the particular test case 
that was being computed. 

5. A temporal preconditioning was added to the 
flow equations to allow solutions at very low 
Mach numbers. The preconditioner was imple- 
mented such that both steady state and time ac- 
curate flows could be computed. Steady state 
solutions were considered in this study. Mach 


number flows as low as 0.0005 were computed 
without degradation to the convergence rate of 
the solution procedure. Without the precondi- 
tioning, convergence was either very slow due to 
the necessity of running at a much smaller time 
step, or the equations could not be converged to 
a solution. Preconditioning was relatively easy 
to add to the numerical code. 
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